clear all;
clc;

tic
disp('Loading data...');
load mnist_test;
load mnist_train;

k = 7;
s = 100;
lam = 1;

C = 10;

X = [xL';xU'];

Y = yL';
Y = [(Y==0) (Y==1) (Y==2) (Y==3) (Y==4) (Y==5) (Y==6) (Y==7) (Y==8) (Y==9)];
Y = Y - (Y==0);
Y = double(Y);
n = size(X,1);
[l,c] = size(Y);
disp('Data loaded!');

toc
tic

disp('Start knn searching...');
[IDX,M] = knnsearch(X,X,'K',k);
disp('knn finished!')

toc
tic

disp('Building graph...');
Ivec = reshape(repmat((1:n)',1,k-1),[],1);
IDXvec = reshape(IDX(:,2:k),[],1);
clear IDX;
Dvec = reshape(M(:,2:k).^2,[],1);
clear M;
Dvec = Dvec/max(Dvec);

A = sparse(Ivec,IDXvec,Dvec);
A = A + A';
A = (A~=0);
disp('Graph finished!');
save('7nngraphplus','A','IDXvec','Dvec');
toc
tic

disp('Check connectivity...');

%check connected
[S,~] = graphconncomp(A,'Directed',false);
if(S > 1) 
    error('Graph not connected! Try more neighbors.');
end
disp('Graph connected!');

toc

B = adj2inc(A);
clear A;

L = double(B)'*double(B);
clear B;
DD = diag(L);
L = 2*sparse(1:n,1:n,DD) - L;
clear DD;

disp('Find eigenvalues...');
tic

[V,~] = eigs(L,s+1,'SM');
save('100eigs','V');
toc

clear L;
V = V(:,1:s);
QC = V(Ivec,:) - V(IDXvec,:);
clear Ivec;
clear IDXvec;

P = zeros(l,s,c);
for i = 1:c
    P(:,:,i) = repmat(Y(:,i),1,s).*V(1:l,:);
end
clear Y;
Q = V(l+1:n,:);
clear V;

disp('Start optimization...');
tic

cvx_begin sdp
cvx_solver sedumi
    variable t
    variable gam(s,c)
    variable M(s,s) symmetric
    variable del(l,1) nonnegative
    minimize (t-lam*trace(M));
    for i=1:c
        [M gam(:,i);gam(:,i)' t-2*C*ones(1,l)*del] >= 0;
        P(:,:,i)*gam(:,i) >= ones(l,1)-del;
    end
    sum(QC*M.*QC,2) <= Dvec
    M == semidefinite(s);
cvx_end

toc

[~,Ypred] = max(Q*gam,[],2);
Ypred = uint8(Ypred)-1;

label = int8(yU');
calacc(Ypred,label)